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Abstract 

Making  use  of  a  500  mb  numerical  forecast  some  problems 
of  atmospheric  large-scale  diffusion  are  investigated.  The 
expansion  of  various  sets  of  clusters  of  particles  is  compared 
with  the  predictions  based  on  the  similarity  theory  of  homogeneous 
turbulence.  The  spreading  rates  were  in  the  beginning  of  the 
period  of  the  forecast  in  a  good  agreement  with  the  theoretical 
predictions,  but  later  on  they  dropped  to  values  considerably 
smaller  than  the  predicted  ones.  The  reasons  of  such  a  behaviour 
are  discussed.  Furthermore  the  difference  in  zonal  and  meridional 
expansions  of  the  clusters  is  demonstrated;  the  zonal  expansion 
rates  are  more  than  one  order  of  magnitude  greater  than  the 
meridional  ones. 

The  intrinsic  property  of  turbulence  to  disturbe  initially 
equidistantly  spaced  particles  until  the  random  distribution 
is  achieved  is  investigated  with  respect  to  its  time  scale. 

It  is  found  that  the  atmosphere  needs  about  6  days  to  move 
the  initially  rectangularly  spaced  particles  381  km  apart 
to  the  positions  which  form  the  distribution  differing 
10  per  cent  from  the  random  distribution. 

1.  Introduction 

It  was  only  comparatively  recently  that  in  diffusion  studies 
the  proper  attention  has  been  paid  to  the  difference  between 
the  diffusion  from  a  fixed  source  and  the  diffusion  relative 
to  a  moving  origin  of  a  single  puff  or  cluster  of  particles. 
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The  former  diffusion  type  was  first  analysed  by  G.I.  Taylor 
in  the  well  known  paper  in  1921*  The  latter  type  was  discussed 
by  L.F.  Richardson  already  in  1926,  but  it  was  not  before 
about  1950  that  in  several  basic  studies  by  Brier  and  Batchelor 
the  detailed  Mathematical  analysis  of  the  problem  was  given* 

Also  a  lot  of  experimental  work  has  been  devoted  to  both 
diffusion  problems*  The  main  achievements  along  all  these 
lines  are  comprehensively  reviewed  in  a  recent  study  by 
Pasquill  (1962)* 

Various  experimental  investigations  cover  the  whole 
range  of  diffusion  processes  which  are  experienced  in  the 
atmosphere,  up  to  the  diffusion  on  the  largest  i.e.  synoptic 
scale*  Here  on  this  synoptic-scale  diffusion  process  our 
attention  will  be  fixed*  A  description  of  the  diffusion  on 
the  synoptic  scale  can  be  achieved  in  a  direct  way  by  the 
observations  of  the  trajectories  of  some  markers  floating 
in  the  air  (constant-level  balloons  for  instance)*  An  alter¬ 
native  procedure  is  to  compute  air  trajectories  in  some  way* 

In  previous  investigations  of  this  sort  (Durst  et  al « ,  1959) 
the  interest  was  concentrated  mainly  in  the  dispersion  of  the 
end-points  of  air  trajectories  initiated  at  regular  intervals  at 
a  fixed  point  (fixed-source  diffusion  type)*  In  this  paper  some 
information  concerning  the  complementary  problem  of  relative 
diffusion  of  clusters  of  particles  will  be  presented. 

The  ideal  way  of  studying  the  expansion  of  a  given 
large-scale  cluster  is  to  know  the  trajectories  of  individual 
particles  of  the  cluster.  To  this  end  the  use  of  constant- 
level  balloons  data  can,  because  of  the  sparsity  of  balloon 
flights  in  space,  hardly  be  made*  Particle  trajectories  can 
be  computed  from  successive  contour  naps  with  the  use  of 
geestrophic,  as  in  the  work  by  Durst  et  al*,'  or  some  other 
wind  approximation*  In  our  study  the  particle  trajectories 
were  obtained  by  the  application  of  a  numerical  weather 
prediction  procedure.  The  main  reason  for  this  choice  lied 
in  computation  facilities,  since  the  numerical  prediction 
program,  suitable  for  the  computation  of  a  large  number  of 


trajectories,  was  available  for  the  investigation  purposes. 

The  use  of  the  **foreca8t'*  instead  of  the  "hindcast"  trajectories 
brings  in  the  diffusion  results  systematic  errors  inherent 
in  the  adopted  forecasting  model  (steady  decrease  in  eddy 
kinetic  energy),  and  this  must  be  borne  in  mind  at  the  exami¬ 
nation  of  the  computation  results. 

The  large  number  of  computed  trajectory  coordinates 
suggested  furthermore  the  idea  to  study  the  time  scale  of 
the  tendency  of  turbulence  to  move  initially  regularly  spaced 
particles  to  the  positions  for  which  the  initial  spacing 
is  completely  disturbed  and  the  random  distribution  is 
achieved.  A  discussion  of  this  problem  is  given  in  section  4. 

As  an  indication  of  this  process  the  mean  value  and  the  fre¬ 
quency  function  of  the  distance  from  every  particle  to  the 
particle  nearest  were  taken. 

2.  Computations  of  the  cluster-expansion  rates 

Trajectory  computations  were  based  on  the  integration 
of  the  set  of  primitive  equations 


VU.CI#.,  -  4^  .  r,  -  I* 


dv  2 . 


o  2  2  2 

2  .  u  -fv  dm: 
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-If 


■(4)r)  , 


valid  for  the  homogeneous  atmosphere  (  Hollmann,  1959  ). 

In  (1)  X  and  y  are  the  Cartesian  coordinates  of  the  image 
surface  of  stereographic  projection,  fr  (  u,v)  is  the  so  called 
^’reduced  image  velocity"  (  Hollmann,  1959  ),  while  other 
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symbols  have  their  conventional  meanings.  Some  details  of 
the  computation  procedure  are  already  described  elsewhere 
(  Mesinger,  1962). 


Fig.  1.  Grid.  Horizontal  grid  distance  381  km  in  the 
image  surface  of  stereographic  projection,  crossing 
the  earth's  surface  at  60^N.  Grid  points  indicated  by 
crosses  are  centers  of  the  clusters  at  the  beginning 
of  the  forecast. 

The  program  available  for  the  solving  of  system  (l) 
on  an  IBM  704  computer  was  designed  for  the  grid  consisting 
of  2080  points  evenly  spaced  on  an  octagonal  area  centered 
on  the  North  Pole.  With  the  adopted  grid  distance  of  381  km 


at  60^N  latitude  the  coaputatioa  area  extended  aonthvard  te 
about  lO^N  latitude  (Fig.  1).  Starting  with  the  geepoteatial 
field  of  the  500  nb  aurfaee  for  the  00  GMT  26  August  1958, 
shown  in  Fig.  2,  a  forecast  for  the  period  of  8  day's  was  nade. 
On'  this  occasion  to  the  kinenatic  eddy-viscosity  coefficient  ^ 
in*(l)  the  value  of  0.517«10  cs  sec  was  prescribed  (roupd 
nujsber  in  binary  systen).  There  is  sene  evidence  that  thlh\ 
couid  be  about  the  right  value  for  the  grid  sise  used  in  the 
cosputations  (Phillips,  1956,  p.l3l). 


Fig.  2.  The  initial  situation:  contour  height  ef  the 
5O0idb  surface,  for  the  00  GMT  26  August  1958. 

Making  use  of  successive  wind  fields  at  intervals  of 
half  an  hour,  obtained  as  an  output  of  this  S-day  forecast. 


the  trajectories  were  coMputed  for  all  the  2080 , particles 
situated  at  the  beginning  of  the  forecasting  period  in  the 
grid  points.  This  was  done  with  the  fomulae  centered  in 
tine  and  space*  Fields  of  particle  coordinates  at  one-henr 
tine  intervals  were  written  on  nagnetic  tapes  as  an  output 
of  the  trajectory  conputation  progran.  Particles  were  then*^ 
conbined  in  suitable  sets  to  fom  various  clusters.  There 
were  52  clusters  of  each  sort,  having  their  centers  at  the 
beginning  of  the  forecasting  period  at  52  .grid  points  evenly 
distributed  over  the  conputation  area  (Fig.  1).  Out  of  the 
tapes  with  the  particle  coordinates  nean  standard  deviations 
of  the  particles  about  the  gravity  centers  of  the  clusters 
(d  'm)  were  evaluated-  at  tine  intervals  of  6  hours.  This 
was  done  with  the  aid  of  two  prograns; 

(a)  A  progran  evaluating  standard  deviations  for  three 
sorts  of  square-fomed  clusters,  consisting  of  3X3,  5X5 
and  7X7  particles,  respectively. 

(b)  A  progran  evaluating  standard  deviations  separately 
in  zonal  and  neridional  directions  for  five  sorts  of  dianond* 
forned  clusters,  consisting  only  of  4  particles  at  1,  2,  3, 

4  and  5  grid  distances  away  fron  the  center,  respectively. 

Since  sone  other  research  on  atnospheric  nacroturbulence 
was  carried  on  parallel  with  these  conputations,  several 
statistical  quantities  useful  for  the  interpretation  of  the 
cluster-expansion  results  were  available  for  this  purpose. 

3.  Discussion  of  the  conputed  cluster-expansion  rates 

Regarding  the  rates  of  expansion  of  the  clusters,  the 
■ost  interesting  point  is  how  far  they  agree  with  the  Batchelor* 
predictions  (  Batchelor,  1950)  based  on  the  siailarity  theory 
of  hoaogeneous  turbulence.  Assuning  that  the  cluster  size 
is  SMll  cospared  to  the  sice  of  the  energy-containing  eddies 
(an  assusption.  which  seems  to  be  fulfiled  at  least  for  the 
SMllest  9-point  clusters)  Batchelor  arrives  at  the 


relations 


^  -  (^^oct^  (small  t),  (2) 

and 

^  ^cC  (intermediate  t),  (3) 

where  the  distinction  between  "small^  and  "intermediate*^ 
times  is  given  approximately  by  the  value 


.2/3  ^-1/3 

%  € 


(4) 


Here  ^  represents  the  average  rate  of  dissipation  of  eddy 

kinetic  energy.  Relation  (3)  should  apply  to  the  clusters 

whose  expansion  is  already  independent  of  the  initial  cluster 

size,  but  is  still  governed  by  the  eddying  motions  belonging 

to  the  inertial  subrange  of  eddy  sizes .  These  formulae  of 

Batchelor  were  tested  by  Gifford  (1957a,  1957b)  on  the  basis 

of  various  sets  of  smoke-puff  spreading  data.  He  concluded 
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that  the  existence  of  the  predicted  t  and  t  regimes  for 
the  growth  of  small  puffs  seems  to  be  strongly  supported  by 
the  experiments. 

To  check  the  applicability  of  the  relations  (2)  and  (3) 

2  2  2 

to  our  large-scale  cluster  data  values  of  ^  -  d*  and  d 

o 

obtained  with  the  program  (a)  were  plotted  vs.  t  on  a  loga- 

2  2  A 

rithmic  paper  (Fig.  3).  Assumption  of  the  law  cx:  t 

gives,  for  the  exponent  oc  values  ranging  for  small  times  from 

1.95  to  2.05;  so  Rel .  (2)  represents  an  excellent  fit  to  the 

2 

data.  An  equivalent  assumption  for  d  gives  for  the  exponent  <K. 
at  the  most  values  ranging  from  1.40  to  1.60.  These  values' 
are  reached  at  times  of  6-8  days,  and  the  inspection  of  the 
trend  of  the  curves  does  not  give  an  indication  that  the 
exponents  could  attain  significantly  greater  values  at  the 
later  times.  Moreover,  exponent  values  should  decrease  after¬ 
wards,  since  in  the  limiting  case  of  very  large  t,  when  the 
whole  energy-spectrum  is  effective,  cluster  must  behave  in 
functionally  the  same  way  as  the  plume,  i.e.  exponent  values 


TIME  (DAYS) 


Fig.  3.  Mean-flquare  particle  diaperaion  (daahed  linea),  and 
increase  of  laean-aquare  particle  diaperaion  relative  to  ita 
initial  value  (full  lines),  for  three  aorta  of  initially 
aquare-fonaed  clusters.  Clusters  consisted  of  9,  29  and  49 
particles,  with  initial  standard  deviations  about  the 
gravity  centers  of  393,  661  and  964  km  respectively. 
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Here  the  questioia  arises  whether  the  8->day  period  is 

long  enough  to  include  the  tine  t^,  after  which  the  t^  regine 

should  be  expected^  To  eTsluate  t^  we  nust  know  the  value 

of  €  f  but  not  nuch  is  known  about  the  right  values  of  g  in 

the  atnosphere*  Since  in  sSope  of  our  investigations  also  a  4-day 

conputation  without  the  energy-dissipation  tern  in  (1)  was  nade^ 

an  estinate  of  e  can  be  nade  fron  the  difference  in  eddy  kinetic 

energies  of  the  two  conputations  (Fig«  4).  During  4  days  the  eddy 

kinetic  energy  (nean  value  over  the  whole  conputation  area) 

of  the  8-day  forecast,  due  to  the  influence  of  the  energy- 

2  -2 

dissipation  tern,  decreased  for  1«32  n  sec  conpared  to  that 
of  the  4-day  forecast*  This  gives  for  the  energy  dissipation 


Fig*  4*  Kinetic  energies  of  the  perturbed  notion  (K*)t  nean 

notion  (K),  and  their  sun,  as  functions  of  tine*  Dashed 

lines  are  the  result  of  the  forecast  nade  without  the  energy- 

dissipation  tern  in  £qa*(l),  and  full  lines  are  the  result 

V  9  2-1 

of  the  forecast  nade  with  the  value  v=  0*517*10  cn  sec 
in  (1). 
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rate  the  value  of  0.038  cm  sec  .  It  is  interesting  that  this 

value  fits  well  in  with  the  trend  in  variation  with  height 

of  several  estimates  of  €  obtained  from  the  small-scale 

diffusion  experiments  (Pasquillf  1962,  p.l61).  Application 

2  “3 

of  the  value  of  0,038  cm  sec  for  €  gives  for  the  time  t^ 
values  ranging  from  4,0  days  (9- point  clusters)  to  7,2  days 
(49-point  clusters);  this  supports  the  idea  that  the  time  t^ 
is  included  in  the  time  scale  of  Fig,  3, 

The  failure  of  Fig.  3  to  show  the  existence  of  the 
t^  regime  should  not  be  considered  as  a  proof  of  inapplica¬ 
bility  of  the  simplif ications  of  inertial  subrange  to  the 
large-scale  atmospheric  diffusion.  We  must  take  into  account 
that  the  forecast  underlying  the  cluster-expansion  computations 
showed  a  steady  decrease  in  eddy  kinetic  energy  (Fig,  4), 
amounted  after  the  8  days  to  41  per  cent  of  the  initial  value. 

This  loss  in  eddy  kinetic  energy  was  due  to  the  (a)  truncation  in 
tine  and  space,  (b)  transformation  of  eddy  kinetic  energy  into  the 
kinetic  energy  of  the  zonally  averaged  motion,  and  (c)  dissipation 
of  eddy  kinetic  energy;  the  dissipation  was  expected  to  differ 
from  its  natural  value  with  respect  to  the  adopted  viscosity 
coefficient.  In  real  atmospheric  stream  fields  the  loss  under 
the  item  (a)  does  not  exist,  and  the  losses  under  the  items 
(b)  and  (c)  are  in  the  average  compensated  by  the  transformation 
of  eddy  potential  energy  into  eddy  kinetic  energy.  Hence, 
in  real  atmospheric  stream  fields  we  must  expect  cluster- 
expansion  rates  to  be  greater  to  a  certain  extent.  This  implies 
that  the  initial  expansion  rates,  found  to  be  in  a  good 
agreement  with  the  Rel . (2) should  also  be  greater. 

Furthermore  we  must  have  in  mind  that  in  our  case  an 
upper  limit  of  the  size  of  the  clusters  exists,  in  contrary 
to  the  considerations  which  led  to  the  prediction  (3),  The 

approaching  to  that  limit  could  also  at  least  partly  be  the 

3 

reason  of  the  lacking  of  the  t  regime* 
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Fig.  5.  Zonallj  averaged  wind,  and  ita  ahaar,  of  tha 
initial  aituation  ahovn  in  Fig«  2« 


In  thaoratical  invaatigationa  of  ralativa  diffuaion 
the  influence  of  the  shear  of  the  aean  flow  is  not  taken 
into  account*  Having  in  nind  the  idea  to  study  this  influence 
experimentally y  in  scope  of  tha  program  (a)  mean  standard 
deviations  were  computed  separately  for  the  40  clusters  south 
of  the  57^N  latitude,  share  the  amount  of  shear  ef  the  mean 
flow  was  comparatively  great,  and  for  the  12  clusters  north 
of  the  57^N  latitude,  where  the  amount  of  shear  was  sbmII 
(Fig.  9).  The  results  of  this  separation  are  shown  in  Fig.  9. 
Clusters  south  of  the  57%  had  greater  spreading  rates  than 
those  north  of  the  97%,  the  difference  being  greater  fer 
clusters  of  larger  initial  sine.  But  since  in  regions  south 
of  the  97%  also  the  eddy  kinetic  energy  was  to  seme  extent 
greater  than  north  of  the  97%  (Fig.  7), it  is  hard  te  say 


i 
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how  much  the  greater  spreading  rates  south  of  the  57  N  were 
a  consequence  of  the  greater  zonal  wind  shear* 


Fig.  6,  Mean  standard  deviations  of  the  40  clusters 
initially  situated  south  of  the  57^N  (full  lines),  and 
of  the  12  clusters  initially  situated  north  of  the 
57*^N  (dot-dashed  lines).  Difference  in  the  initial 
sizes  is  due  to  the  increase  in  grid  distance  with 
latitude* 

In  computations  with  the  program  (b)  only  the  40  clusters 

0 

south  of  the  57  N  were  Included,  since  in  polar  regions  the 
evaluation  of  zonal  and  meridional  distances  with  sufficient 
accuracy  becomes  complicated  and  sometimes  senseless*  The 
number  of  particles  in  the  clusters  was  restricted  to  only 
4  because  of  the  great  consumption  of  time  for  the  automatic 
computation  of  mean  longitude  for  clusters  consisting  of 
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Fig.  7.  Eddy  kinetic  energy  of  the  initial  situation 
shown  in  Fig.  2. 

a  great  nu«ber  of  particles.  Results,  shown  in  Fig.  6,  give 
for  the  zonal  spreading  rates  values  more  than  one  order  of 
magnitude  greater  tlian  for  the  meridional  ones.  It  has  already 
been  described  (Welander,  1955;  Djuric, 1961)  that  the  atmoa* 
pheric  diffusion  has  a  property  to  elongate  fluid  surfaces 
into  long  and  thin  bands.  As  shown  by  this  computation  these 
bands  have  a  strong  tendency  for  an  east-west  orientation. 

We  must  however  mention  that  this  tendency  is  somewhat  exag¬ 
gerated  in  Fig.  8,  since  the  loss  in  eddy  kinetic  energy, 
due  to  deficiencies  of  the  computation  model,  was  much  greater 
for  the  meridional  motion  than  for  the  sonal  one.  The  losses 
in  eddy  kinetic  energies  of  the  meridional • and  the  sonal 
motion  amounted  after  the  8  days  to  66  and  19  per  cent  of 
their  initial  values  respectively. 


! 

i 
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Fig.  8.  Mean  standard  deviations  in  sonal  (full  lines) 
and  in  neridional  (dot-dashed  lines)  direction,  for 
five  sorts  of  initially  diaaond-forsed  clusters. 

The  tendency  for  predoainantly  zonal  stretching  of  the 
clusters  is  further  illustrated  by  the  Fig.  9,  in  which  the 
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Kig.  9.  Shape  of  4  aelectcd  25-point  cluatera  after  the 
8-day  period,  Saae  geonetrical  signs  denote  the  25  par¬ 
ticles  of  a  particular  cluster;  full  sign  of  the  sane 
forn  denotes  their  gravity  center.  Four  squares  represent 
the  initial  positions  of  the  clusters,  and  shaded  areas 
their  approxlnate  8-day  positions.  Full  lines  are  the  f 

trajectories  of  the  gravity  centers  of  the  4  clusters,  | 

and  nunbers  at  their  ends  are  the  8-day  standard  I 

r 

deviations  (in  kilometers)  of  the  particles  of  the  | 

clusters.  The  region  covered  is  the  same  as  in  Figs..  1  and  2. 

form  of  selected  25-poiht  clusters  after  the  8-day  period  is 

demonstrated.  To  make  an  objective  choice  of  the  clust^ers  I 

( 

to  be  demonstrated  in  Fig.  9  all  the  52  25- point  clusters 
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(being  intermediate  to  the  other  two  cluster  sizes)  were 
arranged  according  to  their  8~day  standard  deviation  values, 
and  4  clusters  having  equidistant  positions  in  this  array 
and  at  the  same  time  including  the  extreme  two  clusters  were 
chosen  (  18^^,  35^^  and  52^^  cluster  ).  The  inspection 

of  the  shape  and  orientation  of  the  clusters  shown  in  Fig.  9, 
particularly  that  of  the  uppermost  one  with  the  greatest 
S-day  standard  deviation  value,  shows  that,  at  least  in  our 
computation,  the  upper  limit  to  the  cluster  sizes  is  put 
rather  by  this  zonal  stretching  of  the  clusters  than  by  the 
finite  boundary  of  the  computation  area.  For  a  particular 
cluster  this  limit  is  practically  reached  when  it  completes 
a  zonal  circle,  and  its  center  is  placed  in  the  vicinity  of 
the  pole. 

This  effect  of  approaching  to  the  upper  limit  must 

evidently  influence  mostly  the  greatest  clusters.  Therefore, 

in  order  to  investigate  to  what  extent  this  effect  could  be 

responsible  for  the  lack  of  the  regime,  for  every  sort 

of  square-formed  clusters  three  clusters  with  greatest  8-day 

2 

sizes  were  chosen,  and  their  mean  d  -values  calculated. 

The  results  of  this  calculation  are  shown  in  Table  1. 


Table  1.  Exponent  values  of  the  assumption  O  oc>  t 


Cluster 

type 

3-greatest-clusters  exponents  at 

days 

max.  exponent  i 

1-2 

2-3  3-4 

4-5 

5-6 

6-7 

7-8 

“  considering  all 
the  52  clusters 

7X7 

1.20 

1.70  1.85 

1.80 

1.95 

1.90 

1.90 

1.40 

5X5 

1.25 

2.00  2.30 

1.90 

1.95 

1.80 

2.00 

1.45 

3X3 

1.65 

2.50  2.50 

2.20 

2.25 

2.00 

2.05 

1.60 

The  exponents  of  these  greatest  clusters  of  any  type  reach 
their  maximum  values  somewhat  earlier  than  the  exponents  of 
all  the  52  clusters  together  (Fig.  3).  Later  on  the  exponents  of 
the  greatest  clusters  decrease  a  little,  nevertheless  remaining 
much  greater  than  the  maximum  values  of  the  exponent  of  all 
the  52  clusters.  One  gets  the  impression  that  the  effect  of 
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approaching  the  upper  limit  is  not  significant  for  the  failure 
of  the  computation  to  show  the  t^  spreading  regime. 

4.  On  a  characteristic  tine  scale  of  turbulence 


It  is  the  turbulent  character  of  the  motion  which  gives 
rise  to  continuous  change  of  the  positions  of  moving  particles 
relative  to  each  other.  Thus  when  a  set  of  regularly  spaced 
particles  is  put  into  a  turbulent  notion,  after  a  certain 
time  the  initial  spacing  will  not  more  be  recognizable.  An 
example  of  such  a  situation  is  given  in  Fig.  1  of  the  already 
cited  note  of  the  author  (Mesingcr,  1962),  where  the  2080  end 
points  of  our  S-days  trajectories  are  demonstrated.  In  the 
interior  of  the  computation  area  no  sign  of  regularity  can 
be  seen,  and  the  particles  seem  to  be  distributed  entirely 
at  random.  So  the  questions  arise  (a)  at  what  rate  wilJ  ini¬ 
tially  regularly  spaced  air  particles  under  the  action  of 
atmospheric  motions  approach  to  the  state  of  random  distribu¬ 
tion,  (b)  will  the  atmospheric  motions  realize  such  a  dis¬ 
tribution,  and  (c)  if  they  do  so  how  much  time  do  they  need 
for  that  ?  The  time-scale  of  this  process  is  evidently  depend¬ 
ent  on  the  initial  distance  of  the  particles,  but  for  a  given 
initial  distance  it  is  determined  by  the  intensity  of  turbulent 
motions.  Hence,  this  time-scale  is  an  interesting  character¬ 
istic  of  the  turbulence. 

In  trying  to  answer  these  questions  we  must  first  find 
how  to  determine  to  what  extent  a  given  distribution  of  some 
points  coincides  with  the  random  distribution  ..  A  suitable 
tool  for  this  determination  seems  to  be  the  frequency  function 
of  the  distance  from  every  particle  to  the  particle  nearest 
to  it.  These  **SBallest  distances**,  being  for  the  square  grid 
of  points  all  equal  to  the  grid  size,  must  have  a  definit 
frequency  function  when  the  points  are  distributed  at  random. 
Therefore  our  next  task  will  be  to  derive  the  analytic  ex¬ 
pression  of  this  frequency  function. 
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To  this  end  consider  two  concentric  circles,  and  let 

the  areas  enclosed  by  the  greater  circle  and  by  the  smaller 
2 

one  be  A  and  r  X  respectively.  Now  when  N  points  are  placed 
at  random  in  the  greater  circle,  the  probability  that  no  one 
of  them  will  fall  in  the  smaller  one  is  equal  to  (1  -  r^/A)^, 
The  probability  of  the  complementary  event  that  at  least  one 
point  will  fall  in  the  smaller  circle  is 


1  -  (1  -  r^/A)'^  . 


(5) 


When  we  are  dealing  with  a  great  number  of  points  distributed 
over  a  relatively  large  area  it  is  suitable  to  drop  out  of 
consideration  the  influence  of  the  finite  boundary  of  the 
area  A.  Hence  we  allow  N  and  A  to  tend  to  infinity,  keeping 
at  the  same  time  the  average  number  of  points  per  unit  area 
n  =  N/A  constant.  The  expression  (5)  then  goes  over  to 


and  has  moreover  become  independent  of  the  position  of  the 
center  of  the  ^’smaller*’  circle.  Thus  we  can  imagine  this 
circle  centered  on  any  one  of  the  points  considered,  and  (6) 
will  then  represent  the  probability  P(d  <  r)  that  the  distance  d 
to  the  point  nearest  of  it  takes  a  value  smaller  than  the 
radius  of  the  circle  r.  In  other  words  (6)  is  equal  to  the 
distribution  function  of  d,  and  its  first  derivative 

2 

f(r)  =  2n7rre”^^^  (7) 

is  the  desired  frequency  function  of  the  random  variable  d. 

It  is  also  of  interest  to  know  the  mean  value  of  the 
variable  d,  and  the  mode  of  its  distribution.  Mean  value  is 
equal  to  the  first  moment  of  the  distribution, 

J  =  ^  rf(r)  dr  ,  (8) 

o 

d  =  l/2n.^^^. 


or  according  to  (7),. 


(9) 
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Mode  of  the  distribution  of  d  is  equal  to 

r  =  (2.tCK  (10) 

M 

For  this  value  of  r  the  frequency  function  (7)  takes  its 
■aximuB  value 


f(r  )  =  (2nc>r/e)^/^.  (11) 

ni 

For  a  given  average  density  of  points  it  is  convenient 
to  choose  the  unit  of  length  so  to  have  n  =  1.  Then  we  have 

from  (9),  (10)  and  (ll)  d  =  0.5,  r  =  0.3989  and  f(r  )  =  1.520, 

m  n 

respectively •  Thus,  when  points  making  initially  a  square 
grid  are  brought  to  the  random  distribution,  the  mean  and 
the  most  frequent  distances  from  every  point  to  the  point 
nearest  to  it  amount  to  0.5  and  0,3989  grid  sizes , respect ively . 

On  the  basis  of  the  forecast  described  in  section  2  a 
computation  was  made  with  the  purpose  to  investigate  the 
process  of  approaching  to  the  now  described  random  state. 

The  values  of  d  at  intervals  of  12  hours  were  computed  for 
all  the  970  particles  situated  at  the  beginning  of  the  fore¬ 
cast  in  grid  points  of  the  interior  of  the  left  half  of  the 
grid.  On  this  occasion  the  change  in  grid  distance,  due  to 
change  in  latitude,  was  neglected.  The  choice  of  the  particle 
nearest  to  a  particular  one  ouf  of  these  970  particles  was 
of  course  made  taking  into  account  distances  to  all  the  2079 
remaining  particles. 

Mean  values  of  so  computed  distances  are  shown  in  Fig.  lOr' 
They  can  be  taken  as  a  .rough  indicator  of  the  randomness  of 
the  distribution.  Thus  we  come  to  the  conclusion  that  the 
atmospheric  motions  need  about  6  days  to  realize  a  distribution 
differing  10  per  cent  from  the  random  one,  when  the  initial 
distance  between  particles  is  equal  to  381  km.  The  distribution 
seems  to  approach  asymptotically  to  the  random  distribution. 

The  rate  of  this  approaching  is  less  than  the  exponential 
one;  this  can  readily  be  seen  by  trying  to  find  an  exponential 
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fit  to  the  curve  in  Fig*  10.  Better  fit  to  the  curve  in  Fig,  10 
can  be  realized  by  the  function  of  the  forn 

d  =  0.5  +  0.5(1  +  t/a)"**  ,  (12) 

where  the  appropriate  value  of  b  is  about  2, 


f’i*.  10,  Mean  value  of  distances  fron  every  particle 
initially  situated  in  a  grid  point  to  the  other  such 
particle  nearest  to  it^-as  function  of  time.  The  unit 
of  length  chosen  is  the  grid  size  of  Fig,  1, 

Histograms  of  the  computed  distances  at  time  intervals 
of  2  days I  and  the  frequency  curve  of  the  random  distribution^ 
given  by  £q,  (7),  are  shown  in  Fig*  11*  The  agreement  pt 
computed  histogram  with  the  shape  of  (7)  is  fairly  good  after 
the  6  days,  especially  on  the  decreasing  side  of  the  frequency 
curve.  Computed  height  of  the  first  class  interval  is  at  that 
time  greater  than  the  theoretical  one  after  (7)  because  of 
the  accumulation  of  particles  at  ths  ends  of  the  rows  nearest 
to  the  boundaries  of  the  area  (Mesinger,  1962),  The  random 
state  would  be  approached  more  quickly  without  the  loss  of 


FREQUENCY  FREaUCNCY 


-  21  - 


eddy  kinetic  energy  by  the  nechanisM  described  in  the  preceding 
section » 


Kig.  11.  Uistogrsas  ef  distances  froei  each  particle 
initially  situated  in  a  grid  point  to  the  other  such 
particle  nearest  to  it  (full  lines),  and  frequency 
curves  of  those  distances  for  the  case  of  the  randon 
distribution  of  particles  (dotted  lines).  The  unit  of 
length  chosen  is  the  grid  sise  of  Fig,  1, 

Acknovledgenents,  -  The  conputations  involved  in  sections 
on  cluster  expansion  rates  have  been  perforaed  partly  by 
Miss  0.  Milovanovic,  froa  the  University  of  Belgrade, 

Departaent  of  Meteorology,  The  data  concerning  the  sonal 
profiles  of  the  aean  wind  and  the  eddy  kinetic  energy  (Pigs,  5 
and  7)  were  evaluated  and  kindly  aade  available  to  the  author 
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by  Dr.  P.  Gburcik.  The  program  for  solving  the  system  (1) 
has  been  elaborated  by  the  Research  Section  of  the  Deutscher 
Wetterdienst  in  Offenbach/M* 
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REMARKS  ON  THE  BEHAVIOUR  OF  COMPUTED  AIR  TRAJECTORIES 
NEAR  THE  BOUNDARIES  OF  AN  OCTAGONAL  AREA 


F. 


Mesinger 


In  the  course  of  some  investigations  on  atmospheric 
macroturbulence  carried  out  in  the  Meteorological  Institute 
of  the  Technische  Hochschule  in  Darmstadt,  air  trajectories 
were  computed  for  a  period  of  8  days*  This  computation  revealed 
an  interesting  feature  in  the  behaviour  of  the  trajectories 
near  the  boundaries  of  the  area*  End  points  of  the  trajectories 
showed  a  tendency  to  accumulate  in  the  corners  of  the  computation 
area , 

Trajectory  computations  were  based  on  the  integration 
of  the  primitive  equations  for  the  homogeneous  atmosphere. 
Together  with  the  continuity  equation  they  make  a  complete 
set  of  equations,  which  in  Cartesian  coordinates  of  the  distorted 
plane  of  stereographic  projection  fl]  can  be  written  in  the  form 


3u 

at 


=  m  (  « 


•  7u  +  0  C7  u  )-  - 


2  2  2 

j,  4-  V  dm 

2  Sx^ 


fv 


ax 


=  ro^(  -  + 


+  V  Bm  ^  a<t 

—  ST  ■ 


(1) 


w?S7  •  (4>y) 


In  (l)  v(u,v)  is  the  "reduced  image  velocity"  flj ,  while 
other  symbols  have  their  usual  meaning.  The  system  (1)  was 
solved  with  the  aid  of  a  program  for  the  computer  IBM  704 

^^On  leave  from  the  University  of  Beograd,  Meteorological  Department 
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elaborated  by  the  Research  Section  of  the  German  Weather 
Service.  As  the  boundary  conditions  the  values  ^  (t)  s  O 
at  all  the  boundary  grid  points  were  prescribed.  The  initial 
wind  field  was  obtained  from  the  field  of  geopotentisl  by 
solving  the  balance  equation.  The  forecast  covered  an  octagonal 
area,  centered  on  the  North  Pole  and  extended  southward  to 
about  lO^N  latitude. 

The  program  for  the  solution  of  system  (1)  had  the 

wind  components  in  all  grid  points  and  at  desired  time  intervals 

as  an  output.  Trajectories  for  a  period  of  6  days  were  computed 

for  all  the  2080  grid  points,  with  the  use  of  those  wind 

components  at  time  intervals  of  30  minutes.  The  computation 

was  made  with  centered  time  and  space  differences.  To  the 

9  2  - 1 

viscosity  coefficient  ^  in  (1)  the  value  of  0.517*10  cm  sec 
was  prescribed.  The  computation  was  started  with  the  geopotential 
field  of  the  500  mb  surface  for  the  26  August  1958,  00  GMT. 

An  examination  of  the  end  points  of  the  8-days  trajectories 
(Fig.  1)  showed  an  interesting  feature  concerning  the  spacing 
of  points  in  some  corners  of  the  area.  Let  us  consider  the 
four  rows  of  grid  points  adjacent  to  the  boundaries  of  the 
area  and  parallel  to  the  sides  of  the  figure.  At  three  out 
of  these  four  rows  the  end  points  of  the  trajectories,  which 
started  in  the  grid  points,  moved  westward  along  the  rows 
and  accumulated  at  their  ends.  Moreover,  these  rows  acted 
as  traps  for  the  points  which  reached  them  from  the  interior 
of  the  area.  Fairly  large  regions  without  points  appeared 
due  to  this  unability  of  points  to  leave  those  rows. 

Such  behaviour  of  the  trajectories  appeared  as  a  conse** 
quence  of  the  fact  that  at  these  four  rows  the  velocity  com¬ 
ponents  normal  to  the  boundary  were  extremely  small  during 
the  whole  period  of  the  forecast.  In  the  beginning  these 
normal  velocity  components  were  equal  to  zero:  since  the 
initial  wind  field. was  obtained  by  the  use  of  the  stream 
function,  and  both  velocity  components  had  to  be  zero  at 
the  boundary  points,  a  constant  value  had  to  be  prescribed 
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Fig.1.  End  points  of  th«  8-days  trajoctorias  which  originotsd 

In  the  grid  points 

to  the  stream  function  at  the  boundary  points  and  the  points 
adjacent  to  the  boundary.  As  a  consequence  the  velocity  coapo 


if* 
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nents  normal'  to  the  boundary  were,  also  zero  at  the  four  con¬ 
sidered  rows  adjacent  to  the  boundary.  During  the  first  two 
days  of  the  forecast  these  noraal  velocity  components  increased 
only  to  about  0,2  m  sec  after  that  time  they  remained  at 
this  small  order  of  magnitude. 

The  accumulation  and  depletion  of  the  trajectory  points 
is  likely  to  be  accompanied  with  the  mass  accumulation  and 
depletion.  Numerical  instability  which  appears  always  after 
a  certain  time  at  the  boundaries  of  the  computation  area 
could  be  a  consequence  of  this  mass  accumulation  and  depletion. 

Acknowledgement:  The  program  for  solving  the  system  (1) 
has  been  elaborated  by  the  Research  Section  of  the  Deutscher 
Wetterdienst  in  Offenbach/M. 
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